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ABSTRACT 

A  strain-based  thermodynamics  framework  is  proposed  for  modeling  the  continuum 
damage  behavior  of  viscoelastic  materials.  Damage  is  represented  by  an  internal 
state  variable  in  the  form  of  a  symmetric  second  rank  tensor.  The  effect  of  damage  on 
the  constitutive  behavior  is  introduced  through  direct  coupling  between  the  damage 
variable  and  the  viscoelastic  internal  state  variables.  This  approach  accounts  for 
time-dependent  damage  as  well  as  damage-induced  changes  in  material  symmetry. 
Also,  damage  evolution  is  modeled  by  employing  the  concept  of  damage  surfaces. 


1 


1.  INTRODUCTION 


The  growing  interest  in  the  use  of  polymeric  materials  (e.g.  plastics  and  polymeric 
composites)  for  structural  applications  mandates  appropriate  knowledge  of  the  me¬ 
chanical  behavior  as  well  as  the  durability  of  these  materials.  It  is  well-known  that 
polymeric  materials  creep  viscoelastically.  In  addition,  experimental  investigations 
(e.g.  Suvorova,  1985;  Schapery,  1989;  Tuttle  et  a/.,  1995;  Smith  and  Weitsman,  1996) 
indicate  that  polymeric  composites  may  undergo  distributed  damage  in  the  form  of 
a  multitude  of  microcracks.  The  objective  of  this  article  is  to  establish  a  framework 
for  the  constitutive  modeling  of  the  foregoing  features  of  material  behavior.  Such 
framework  is  essential  for  a  reliable  engineering  design. 

Up  to  the  present  time,  most  of  the  efforts  for  modeling  distributed  damage 
have  been  directed  toward  brittle  materials  exhibiting  elastic  behavior  and  metals 
exhibiting  plastic  or  creep  response  (e.g.  Krajcinovic,  1996).  Less  attention  has 
been  paid  to  the  modeling  of  damage  in  viscoelastic  materials.  Notably,  Schapery 
(1981,  1994,  1996a)  established  a  basic  formulation  for  viscoelastic  response  that 
is  accompanied  by  microstructural  changes,  such  as  profuse  micro  cracking.  The 
microstructural  changes  are  represented  in  Schapery’s  work  by  means  of  a  set  of 
internal  state  variables  whose  evolutionary  laws  are  motivated  by  considerations  of 
viscoelastic  fracture  mechanics.  It  may  also  be  mentioned  that  Weitsman  (1988) 
attempted  to  model  the  coupling  between  viscoelasticity  and  damage  for  a  special 
class  of  linear  viscoelastic  materials. 

In  a  recent  article  by  the  present  authors  (Abdel-Tawab  and  Weitsman,  1996) 
a  stress-based  formulation  for  modeling  the  coupling  between  viscoelasticity  and 
distributed  damage  was  developed  and  applied  to  a  swirl-mat  polymeric  composite. 
The  effect  of  damage  on  material  behavior  was  introduced  through  the  concept  of 
effective  stress.  Also,  damage  evolution  was  related  by  the  empirical  Kachanov- 
Rabotnov  forms  (Kachanov,  1986),  which  are  best  suited  for  monotonic  creep  load¬ 
ings.  For  more  complex  loading  histories  the  concept  of  damage  surfaces  (Krajci¬ 
novic,  1996)  offers  a  more  versatile  approach  to  damage  evolution. 
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It  is  well-known  that  damage  surfaces  are  better  expressed  in  strain  space  than 
in  stress  space  (Ju,  1989;  Krajcinovic,  1996).  This  concept,  which  seems  particularly 
appropriate  for  viscoelastic  response  -  where  creep  occurs  at  all  stress  levels,  pro¬ 
vided  a  motivation  for  the  present  strain-based  formulation  of  viscoelasticity  coupled 
with  damage.  An  additional  motivation  is  the  fact  that  strain-based  viscoelastic  con¬ 
stitutive  models  are  more  convenient  for  implementation  into  finite  element  codes 
than  stress-based  ones  (e.g.  ABAQUS,  1996).  The  present  formulation  employs 
concepts  of  continuum  damage  mechanics  as  well  as  several  existing  concepts  of  the 
thermodynamic  theory  of  viscoelastic  materials  (Biot,  1954;  Schapery,  1964).  This 
format  accounts  for  time-dependent  damage  as  well  as  damage  induced  changes  in 
material  symmetry. 

In  Section  2  of  this  article,  we  present  a  general  thermodynamics  framework  that 
accounts  for  both  viscoelastic  and  damage  processes.  We  proceed  by  modeling  the 
coupling  between  these  two  processes  in  Section  3.  In  Section  4,  damage  evolution 
is  modeled  through  the  concept  of  damage  surfaces  and  is  illustrated  by  a  simple 
example  in  Section  5.  Section  6  concludes  with  a  summary  and  some  remarks 
pertinent  to  the  present  work. 

2.  THERMODYNAMICS  FRAMEWORK 

Consider  a  polymeric  material  and  let  7^  (r  =  1,2,. . .  ,R)  denote  R  scalar  valued 
internal  state  variables  representing  the  internal  degrees  of  freedom  of  molecular  mo¬ 
tion  in  the  polymeric  chains.  The  internal  state  variable  representing  damage  can 
be  related  in  terms  of  tensorial  quantities  of  even  ranks,  which  can  be  associated 
with  the  spatial  distributions  of  microcracks  (Krajcinovic,  1996).  For  simplicity,  the 
damage  variable  is  chosen  as  a  symmetric  second  rank  tensor  a;, 7  with  dimensionless 
components.  This  damage  variable  is  capable  of  simulating  changes  in  material  sym¬ 
metry  such  that  an  initially  damage  free  isotropic  material  may  become,  at  most, 
orthotropic  upon  damage  formation  (Cordebois  and  Sidoroff,  1982;  Lemaitre,  1992). 
Despite  the  shortcomings  of  the  abovementioned  damage  variable  (e.g.  Lubarda  and 
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Krajcinovic,  1993),  it  was  adopted  by  several  workers  in  the  field  of  damage  me¬ 
chanics  (e.g.  Cordebois  and  Sidoroff,  1982;  Murakami  and  Imaizumi,  1982;  Chen 
and  Chow,  1995)  due  to  its  relative  simplicity  and  applicability  to  practical  circum¬ 
stances.  It  should  be  mentioned  that  the  present  formulation  can  be  readily  modified 
to  accommodate  damage  variables  of  other  tensorial  ranks.  Throughout  this  arti¬ 
cle  the  subscripts  r  and  q  are  reserved  for  scalar  quantities;  and  a,  b,  c,  d,  i,j,  k,  I,  m 
and  n  are  associated  with  tensorial  quantities  and  cover  the  range  1,2,3.  Also,  the 
summation  convention  is  implied  over  the  range  of  repeated  indices  unless  stated 
otherwise. 

Viscoelasticity  and  damage  are  irreversible  thermodynamic  processes.  For  a 
closed  system  and  small  strains,  the  entropy  production  inequality  can  be  written 
in  the  form  (Coleman  and  Curtin,  1967) 

+  >0,  (1) 

where  xj)  is  the  Helmholtz  free  energy  (per  unit  volume),  <Jij  -  components  of  a 
suitably  defined  volume  averaged  stress  tensor,  Cij  -  components  of  the  infinitesimal 
strain  tensor,  S  -  entropy  (per  unit  volume),  T  -  temperature,  hi  -  components  of 
the  heat  flux  vector,  T,  =  dTjdxi  -  components  of  the  temperature  gradient,  and 
Xi  -  space  coordinates.  Also,  in  (1)  an  overdot  signifies  differentiation  with  respect 
to  time. 

Consider  a  Helmholtz  free  energy  of  the  form 


xp  =  xl>{eij,'yr,iOab,T)  . 


(2) 


The  function  xp  is  assumed  to  be  continuous  and  sufficiently  differentiable  with 
respect  to  its  arguments.  Considerations  of  the  entropy  production  inequality  in  (1) 
together  with  the  functional  dependence  in  (2)  give  the  following  familiar  relations 


a 


ij 


dxp 
dsij  ’ 


(3) 


<S  = 


dxp 


4 


(4) 


and 


h  T 

Trir  +  ^ab<^ab-^  >  0, 


(5) 


where  and  flab  are  the  thermodynamic  forces  conjugate  to  the  internal  state 
variables  7^  and  uJab,  respectively,  and  are  given  by 

dtp 


and 


r.  = 


57.’ 

dtp 

dUab 


(6) 


(7) 


Finally,  from  the  dissipation  inequality  (5)  we  have  the  following  requirements 


r.7r>0,  (8) 

Fr  7r  +  ^  0.  (9) 

Inequality  (8)  must  be  satisfied  whenever  viscoelastic  deformation  occurs,  while 
when  deformation  is  accompanied  by  damage  inequality  (9)  should  be  satisfied  as 
well. 


3.  CONSTITUTIVE  MODELING 

3.1.  General  Formulation 

In  this  section  a  constitutive  model  is  formulated  for  the  case  of  linear  viscoelastic 
behavior  coupled  with  damage.  For  simplicity,  attention  is  restricted  to  the  case  of 
isothermal  behavior.  The  extension  to  the  general  case  of  nonisothermal  conditions 
can  be  made  following  the  same  approach  adopted  here.  The  formulation  will  be 
first  established  for  fixed  strain  Sij  and  damage  (jJab  and  subsequently  extended  to 
fluctuating  Cij  and  ujab- 

For  fixed  Sij  and  Uab  an  irreversible  thermodynamic  process  is  triggered  in  the 
material,  which  prompts  the  viscoelastic  internal  state  variables  jr  to  drift  sponta¬ 
neously  toward  their  equilibrium  values  7®.  Under  isothermal  conditions,  all  7®  are 
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independent  of  temperature,  hence 


(10) 


These  equilibrium  values  are  assumed  to  be  continuous  and  sufficiently  differentiable 
functions  of  their  arguments.  Assuming  that  all  7^  and  7®  are  sufficiently  small,  a 
Taylor  series  expansion  for  ^  about  7®  takes  the  form 


where 


V’  =  V’e  +  ^  V’r?  (7r  “  7r)  {iq  “  7«)  +  H.O.T.  , 


(11) 


is  the  value  of  V’  at  equilibrium, 


'^rq 


\d^rd^q  )  g 

is  a  symmetric  matrix  considered  to  be  constant,  and  H.O.T.  refers  to  higher  order 
terms  neglected  due  to  smallness  of  7^  and  7®.  In  the  above  relations,  and  in  the 
sequel,  the  subscript  “  e  ”  implies  that  a  quantity  is  calculated  at  7r  =  7®  V  r.  Note 
that  at  equilibrium  ‘tp  is  minimum  (Callen,  1960;  Prigogine,  1967),  and  hence 


and 

'Pvq  ^  0  . 

Consequently,  there  is  no  linear  term  in  (11)  and  xprq  is  a  positive  definite  matrix.  It 
should  be  mentioned  that  an  expansion  similar  to  that  in  (11)  was  previously  used 
by  Lubliner  (1972). 

Employing  the  usual  assumption  of  viscous-like  resistance  (Biot,  1954;  Schapery, 
1964),  let 

Tr  ~  ®rg  7g  ?  (f 
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where,  according  to  Onsager’s  principle  (Callen,  1960;  Fung,  1965),  is  a  sym¬ 
metric  matrix.  Substitution  of  (12)  into  inequality  (8)  gives 

^rq  'jr  'jq  ^  0  . 

Hence,  the  matrix  is  positive  semi-definite.  Note  that  in  the  general  case  a^q  is 
a  function  of  temperature,  but  since  we  are  considering  only  isothermal  conditions 
then  Crg  is  constant. 

Equations  (6),  (11)  and  (12)  yield 


®rg  ifq  0rg  'Jq  V’rg  'Jq  ■  (1-^) 

Since  a^q  is  a  constant  symmetric  positive  semi-definite  matrix  and  if^rq  is  a  constant 
symmetric  positive  definite  matrix,  it  is  possible  to  rewrite  (13)  in  a  diagonalized 
form  (Meirovitch,1967)  as 

7r  =  It  (>“0  sum  over  r)  ,  (14) 

where  7r  are  transformed  internal  state  variables,  each  being  a  linear  combination  of 
the  original  internal  state  variables  7^.  The  parameters  7®  are  the  equilibrium  values 
corresponding  to  7^  and  are  obtained  from  7®  by  the  same  linear  transformation  as 
that  for  7r.  Also,  At  and  are  constants  such  that  >  0  and  >  0. 

For  fixed  strain  and  damage,  the  solution  of  equation  (14)  is 

7r  =  7r  (1  —  (no  sum  over  r) ,  (15) 

where  are  relaxation  times  given  by 

Tt  =  (no  sum  over  r)  .  (16) 

In  terms  of  the  transformed  internal  state  variables,  expansion  (11)  can  be  rewritten 
as 

=  +  i',f  +  H.O.T.  (17) 
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The  viscoelastic  strain  can  now  be  obtained  by  substituting  (17)  into  (3)  bearing  in 
mind  that  7^,  and  hence  7^,  are  to  be  kept  fixed  during  the  partial  differentiation 
indicated  in  (3).  Employing  (15)  we  then  obtain 


dfpe  .  dAj-  ^-ifTT 
deij  r  ^ 


(18) 


where 


Ar  =  Ar{eij,U!ab)  =  \'^r  (7r)^ 


(no  sum  over  r) 


(19) 


The  first  term  on  the  right  hand  side  of  (18)  represents  the  long-term  (rubbery) 
part  of  the  behavior,  and  the  second  term  represents  the  transient  (time-dependent) 
part. 

Motivated  by  previous  works  on  linear  elasticity  with  damage  (e.g.  Lemaitre 
and  Chaboche,  1985;  Lemaitre,  1992),  we  now  recast  the  formulation  in  a  format 
that  retains  a  linear  viscoelastic  relaxation  modulus  and  introduces  the  effects  of 
damage  by  mapping  the  stress  and  strain  into  “damage  effective”  stress  and  strain, 
respectively.  To  this  end,  consider  first  the  transient  part  in  (18).  Expanding  Ar  in 
terms  of  strain  around  the  reference  state,  =  0,  up  to  quadratic  terms  to  retain 
linearity  one  obtains 


A  _  1 

Ar  -  2 


d'^Ar 


^ij  ^kl  1 


(20) 


/  0 

where  the  subscript  0  implies  that  a  quantity  is  calculated  at  the  reference  state. 
Note  that  the  constant  term  in  (20)  vanishes  since  (7®)o  =  0,  and  hence  7®  =  Ar  =  0 
at  the  reference  state.  In  addition,  the  linear  term  in  (20)  is  discarded  since  it 
corresponds  to  a  residual  stress  at  the  reference  state,  which  is  disregarded  in  the 
present  formulation. 

A  more  specific  functional  form  for  Ar  can  be  obtained  by  realizing  that  the 
internal  molecular  motions  represented  by  7r  occur  on  a  much  smaller  dimensional 
scale  than  that  of  damage  represented  by  ujab-  This  suggests  that  all  7®,  and  hence 
all  7®  and  Ar,  are  likely  to  be  affected  by  damage  in  a  common  manner;  i.e.  they 
have  common  dependence  on  Uab-  Consequently,  we  can  rewrite  (20)  in  the  form 


Ar  —  2  ^ija-b  abed  Pcdkl  ^ij  ^kl  V  T  . 


(21) 
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In  (21),  is  a  double  symmetric  fourth  rank  tensor  (i.e.  = 

ACljik  =  ACliij),  and  Pijki  =  Pijkiii^ab)  is  a  double  symmetric  fourth  rank  tensor 
valued  function  of  the  damage  variable  ujab  such  that 

at  (-^ab  ~  0  ^  P ijkl  ~  Pjkl  ? 


where 

^ijkl  2  i.^ik  ^jl  d”  Sil  Sjk^  , 

is  the  unit  fourth  rank  tensor  and  Sij  is  Kronecker  delta. 

Using  (21),  the  transient  part  of  (18)  takes  the  form 

^  =  Pijab  ACabcd{t)  Pcdkl  &kl  ,  (23) 

r  OSij 


where 

=  Y.  AC'h  e-""  .  (24) 

r 

From  (18),  (22)  and  (23)  it  is  clear  that  ACiju  is  the  undamaged  transient  (time- 
dependent)  stiffness  tensor. 

Consider  now  the  long-term  part  of  (18).  The  equilibrium  Helmholtz  free  energy 
V’e  can  be  expanded  around  the  reference  state  in  the  form 

,  1  (  \ 

’*•“2 

Following  common  practice  in  continuum  damage  mechanics  (e.g.  Lemaitre  and 
Chaboche  1985;  Lemaitre,  1992),  we  further  assume  that  depends  on  damage  in 
the  separable  form 

V’e  —  2  ^abcd  Qcdkl  ^ij  ^kl  t  (26) 


where  Qijki  =  Qijkii'^ab)  is  a  double  symmetric  fourth  rank  tensor  valued  function 
of  u^ab  such  that  Qijki  =  hjki  at  Uab  =  0,  and  consequently  Clji.i  is  the  undamaged 
long-term  (rubbery)  stiffness  tensor  which  is  typically  positive  definite  (Fung,  1965). 
Using  (26),  the  long-term  part  of  (18)  takes  the  form 


dsij 


Qijab  (^abcd  Qcdkl  ^kl  ? 


(27) 
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Relations  (23)  and  (27)  indicate  that  the  long-term  part  of  the  behavior  can  in 
general  depend  on  damage  in  a  manner  that  differs  from  that  of  the  transient  part. 
However,  for  simplicity,  we  assume  here  that  both  parts  have  the  same  dependence 
on  damage  so  that 

Qijki  =  Pijki  ■  (28) 


Assuming  that  the  inverse  P-h  exists,  define  the  following  “damage  effective”  stress 
and  strain  tensors 


dij  =  -  aki  . 

(29) 

Sij  —  Pijkl  ^kl  ? 

(30) 

then  relation  (18) 

can  be  rewritten  in  the  compact  form 

^ij  ^ijklij'^  ^kl  ? 

(31) 

where 

Cijki{t)  =  +  ACijki{t) , 

(32) 

is  the  overall  (long-term  and  transient)  stiffness  tensor  of  the  undamaged  material. 
Relations  (29)  and  (30)  are  consistent  with  the  formulations  of  the  concepts  of 
effective  stress  and  effective  strain  (e.g.  Cordebois  and  Sidoroff,  1982;  Simo  and  Ju, 
1987),  where  the  mapping  tensor  for  the  effective  stress  is  taken  to  be  the  inverse  of 
that  of  the  effective  strain.  Also,  note  that  both  aij  and  are  symmetric  due  to 
the  hypothesized  double-symmetry  of  Pijki- 

Relation  (31)  suggests  that  for  a  given  damage  level  Uabi  instantaneous  mapping 
of  the  actual  stress  (Tki  and  strain  Ski  according  to  (29)  and  (30),  respectively,  lead  to 
new  stress  and  strain  iij  quantities  that  are  related  by  the  usual  linear  viscoelastic 
constitutive  relation  for  fixed  strain  (e.g.  Fung,  1965).  Upon  hypothesizing  time- 
translation  invariance,  and  since  aij  is  linear  in  iki-,  a  straightforward  application  of 
the  superposition  principle  (Pipkin,  1986)  to  expression  (31)  yields 

=  f  Cijki{t  -  t)  dr  .  (33) 

Jo-  cLt 
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Allowing  for  spatial  variations  of  stress  and  damage  the  total  derivative  dj dr  inside 
the  integral  is  replaced  by  a  partial  derivative  dj 5r,  holding  the  spatial  coordinates 
Xi  fixed.  Thus 

dij  =  J^_  Cijki{t  -  t)  •  (34) 

Equation  (34)  is  the  stress-strain  constitutive  relation  for  the  coupled  linear  vis¬ 
coelastic/damage  behavior,  and  can  be  expressed  in  terms  of  the  actual  stress  and 
strain  as 

(TiJ  =  Pijab  jl_  Cabcdit  -  t)  .  (35) 

It  should  be  noted  that  at  t  =  0 

c««(0)  =  CJh  =  CJh  +  ,  (36) 

r 

where  is  the  initial  (elastic)  stiffness  tensor  which  is  positive  definite  (Fung, 
1965).  For  the  special  case  of  isotropic  virgin  material  response,  the  overall  stiffness 
tensor  Cijki  takes  the  form  (Fung,  1965) 

Ci^kiit)  =  2  G{t)  h^ki  +  [m  -  ^  C?(f)]  Sij  Ski ,  (37) 

where  G{t)  is  the  overall  shear  modulus  and  K{t)  is  the  overall  bulk  modulus  given, 
respectively,  by 

(?(t)  =  Go  +  AG(t) ,  (38) 

and 

K{t)  =  K,  +  A/f (t) .  (39) 

In  the  above  expressions  Go  and  Ko  are  the  instantaneous  shear  and  bulk  mod¬ 
uli,  respectively;  and  AG{t)  and  AK{t)  are  the  transient  shear  and  bulk  moduli, 
respectively,  obtained  from  (24)  as 


AG(t)  =  , 

r 

(40) 

AK{t)  =  , 

(41) 

r 
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where  AGr  and  AKr  are  positive  constants. 


3.2.  The  Dissipation  Inequality 

The  thermodynamic  force  Clab  conjugate  to  ooab  can  be  obtained  by  substituting 
(17)  into  (7)  using  (21),  (26)  and  (28) 


^  _  p  ^  ^Pmnkl 

—  ^ijcd  ^  cdmn  ^ 


ab 


(42) 


where  Cijki  is  given  by 


=  ct,„  -  E  I  I  . 


7r 


(43) 


Since 


7 


0  <  <  1;  Vr, 

V7r/ 


then  Cijki  is  bounded  by  and  Cfjf.i  corresponding,  respectively,  to  the  upper 
and  lower  limits  of  7r/7r-  Since  Cfji^i  and  Ciji.i  are  positive  definite  then  it  follows 
that  Cijki  is  also  positive  definite. 

Employing  (6)  and  (7),  the  dissipation  inequality  (9)  can  be  expressed  as 

;;2  ^ 

/  -  7r  Pijab  Cabod  P :dkl  ^ij  ^kl  ^  0  i 


where 


•  dPcdki  . 

J  zdkl  —  O  ^7i 


(44) 


(45) 


Noting  that  the  first  term  on  the  left-hand  side  of  (44)  is  always  non-negative,  then 
a  sufficient  but  not  necessary  condition  to  satisfy  (44)  is 


P ijab  Cabcd  P tdkl 


positive  semi  —  definite  . 


(46) 


3.3.  The  Mapping  Tensor 

The  functional  form  of  the  mapping  tensor  Pijki  is  restricted  by  the  requirement 
that  Pijki  is  double  symmetric  in  addition  to  the  requirements  in  (22)  and  (44).  In 
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general,  Pijki  is  an  anisotropic  fourth  rank  tensor  function  of  ujab-  However,  due 
to  the  complexity  of  anisotropic  functional  forms  (Zheng,  1994)  and  the  fact  that 
damage-induced  anisotropy  (or,  more  precisely,  orthotropy)  can  be  deduced  from  a 
symmetric  second  rank  damage  tensor,  then  the  more  complex  anisotropic  functional 
form  may  be  avoided. 

Following  Murakami  and  Imaizumi  (1982),  a  simpler  representation  of  Pijki  can 
be  obtained  by  taking  it  as  an  isotropic  fourth  rank  tensor  function  of  iVab-  A  further 
simplification  of  the  representation  of  Pijki  is  obtained  by  considering  a  case  of  dilute 
concentration  of  microcracks  in  which  Pijki  is  linear  in  uJab-  In  this  case,  Pijki  can 
be  written  as  (Murakami  and  Imaizumi,  1982) 

Pijki  =  Cl  Sij  Ski  +  C2  {Sik  Sji  -h  Sil  Sjk)  +  C3  Sij  LCkl  +  C4  Ski  ‘^ij 

-f  C5  {Sik  ^ji  +  ^ii  ‘^jk  +  Sjk  ojii  +  Sji  ujik)  +  H.O.T. ,  (47) 

where  (a  =  1, 2, . . . ,  5)  are  constants. 

From  (47),  it  is  clear  that  double  symmetry  of  Pijki  dictates  that  C3  =  C4.  Also, 
the  requirement  in  (22)  renders  Ci  =  0  and  C2  =  1/2.  To  determine  C3  =  C4  and 
C5,  we  consider  the  special  case  of  isotropic  damage  in  which  the  ensuing  damage 
pattern  does  not  affect  the  symmetry  of  the  virgin  material.  In  this  case  damage  is 
represented  by  a  single  scalar  u  so  that 

<^ab  =  ^Sab  ;  0  <  W  <  1  ,  (48) 

and  the  mapping  tensor  Pijki  takes  the  form  (Simo  and  Ju,  1987) 

Pijki  =  (1  —  n;)  lijki ,  (49) 

which  is  the  inverse  of  the  corresponding  mapping  tensor  that  maps  the  applied 
stress  into  the  Kachanov  effective  stress  in  the  case  of  scalar  damage  (e.g.  Kachanov, 
1986).  The  functional  form  in  (49)  can  be  recovered  from  (47)  by  setting  C3  =  C4  =  0 
and  taking  C5  =  —1/4.  Thus,  the  simplest  possible  form  of  Pijki  becomes 

Ptjkl  =  Pjkl  ~~  {Sik^jl  +  Sil  LOjk  -|-  Sjk  l^il  -1-  Sji  U>ik)  .  (50) 
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It  is  interesting  to  note  that  P^h  obtained  from  (50)  indeed  coincides  with  one  of 
the  forms  proposed  by  Chen  and  Chow  (1995)  for  the  tensor  that  maps  the  applied 
stress  into  an  effective  stress. 

The  complete  formulation  of  the  constitutive  model  requires  an  expression  for 
the  evolution  of  the  damage  tensor  ioab  such  that  (44)  is  satisfied.  Such  an  expres¬ 
sion  can  be  formally  derived  from  thermodynamic  considerations  (e.g.  Lemaitre  and 
Chaboche,  1985;  Lemaitre,  1992),  but  the  usefulness  of  such  approach  seems  to  be 
restricted  to  elastic  response  with  damage.  In  practice,  the  form  of  the  damage  evo¬ 
lution  equation  depends  on  the  material  considered  and  the  applied  loading.  This 
dependence  is  better  correlated  within  the  concept  of  damage  surfaces  (Krajcinovic, 
1996)  as  discussed  in  the  following  section. 


4.  DAMAGE  EVOLUTION 

The  approach  adopted  here  for  describing  damage  evolution  follows  closely  that  pre¬ 
sented  by  Simo  and  Ju  (1987)  and  Lubarda  and  Krajcinovic  (1995).  This  approach 
has  two  main  ingredients.  First,  a  damage  surface  is  introduced  in  strain  space  to 
distinguish  between  the  material  states  associated  with  evolving  damage  and  those 
with  stationary  damage.  Second,  a  damage  potential  is  assumed  to  exist,  from  which 
the  constitutive  law  of  damage  growth  (i.e.  the  damage  rate  Coab)  can  be  derived. 

To  characterize  damage  evolution,  i.e.  damage  loading  conditions,  a  damage  func¬ 
tion  f{eij,K)  is  introduced  so  that 


/(£»;,«)  <0,  (51) 

where  ac  is  a  positive  scalar  damage  threshold  history  parameter  and  at  the  initial 
onset  of  damage  k  —  Ko-  The  equality  in  (51),  i.e.  /  =  0,  corresponds  to  strain  states 
that  lie  on  the  damage  surface  and  for  which  damage  can  evolve.  For  simplicity,  the 
function  /  is  chosen  in  the  simple  isotropic  hardening  form 

f{e^j,K)  =  P{€ij)  -  K  ,  (52) 
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where  ^  is  a  scalar  function  of  the  strain. 

Introduce  a  monotonic  scalar  function  Q{^ab)  such  that  the  damage  rate  can  be 
expressed  as 

<^a6  =  ^ 

Ollab 

where  A  is  a  monotonically  increasing  positive  scalar,  i.e. 

A  >  0  .  (54) 

Physically,  A  represents  a  measure  of  the  cumulative  damage  at  the  considered 
instant  of  the  deformation  process.  The  function  Q  is  referred  to  as  the  “damage 
potential” . 

Following  Simo  and  Ju  (1987)  let 

A  =  /c ,  (55) 

and  define  damage  loading/unloading  conditions  according  to  relations  (51)  and  (54) 
together  with 

A/  =  0.  (56) 

Thus,  if  /  <  0  then  A  =  0  and  from  (53)  no  damage  evolution  takes  place,  i.e.  the 
so-called  damage  unloading  from  the  current  state  of  strain  on  the  damage  surface 
takes  place.  If  /  =  0  and  A  =  0  then  damage  neutral  loading  occurs.  Finally,  if 
A  ^  0  then  /  =  0  and  damage  loading  takes  place. 

During  damage  loading,  the  consistency  condition 

/  =  0  ,  (57) 

must  always  be  satisfied.  From  (52)  and  (57)  we  have 

.  ■ 

=  (58) 

Assuming  that  no  damage  healing  occurs,  i.e.  the  damage  surface  can  only  expand, 
then  K  is  obtained  from  (52)  and  (55)  as 

K  —  max  "(aCq,  .^inax  }  ^  (b9) 
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where  J^raax.  is  the  maximum  value  of  T  over  the  entire  loading  history.  Substitution 
of  (58)  and  (55)  into  (53)  yields 


.  _  . 

d^ah  deij 


(60) 


Thus,  specification  of  the  functional  forms  of  J-  and  Q  completes  the  formulation  for 
the  damage  evolution.  In  practice,  these  functional  forms  depend  on  the  material 
considered  and  the  ensuing  damage  pattern.  Example  functional  forms  will  be 
presented  in  the  following  section. 

In  the  damage  evolution  equation  (60)  the  thermodynamic  force  Q.ab  is  given 
by  expression  (42)  in  which  Cijki  is  given  by  (43).  Thus  an  explicit  expression  for 
Q.ah  requires  evaluation  of  the  ratio  !  Ir-  This  ratio  can  be  determined  from  the 
differential  equation 


d  ( 1  d  Kr\  ( ^r\  1 

dt  \7'/  \Tr  2Ar  dt  )  \7®/  Tr 


(no  sum  over  r) ,  (61) 


where  Ar  is  given  by  (21).  Equation  (61)  is  obtained  after  simple  algebraic  manip¬ 
ulations  of  equation  (14)  and  making  use  of  (16)  and  (19). 

It  should  be  noted  that  D,cd  depends  on  ujab  explicitly  through  Pmnki  and  also 
implicitly  through  the  ratio  7r  /  7r  •  Thus,  in  practice,  equations  (34),  (60)  and  (61) 
need  be  implemented  incrementally  where  for  given  strain  and  time  increments  an 
iterative  procedure  is  required  for  determining  the  corresponding  damage  increment. 


5.  ILLUSTRATIVE  EXAMPLE  FOR  THE  DAMAGE  FUNCTIONS 


Two  scalar  functions  P'{eij)  and  G{0,ab)  are  needed  in  (60)  to  obtain  an  explicit 
damage  evolution  relation.  The  simplest  possible  representation  of  these  functions 
is  to  take  each  as  an  isotropic  function  of  its  tensor  argument.  Thus 


(62) 

(63) 
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where  /f,/2  isotropic  invariants  of  the  corresponding  second  rank 

symmetric  tensor  ^ij.  These  isotropic  invariants  can  be  written  as  (Zheng,  1994) 

4  =  0kk , 

and 

/f  =  det[Aj] , 

where  I3\j  is  the  deviatoric  part  of  /Sij 

To  simplify  matters,  discard  the  dependence  of  ^  on  /|  -  this  is  a  customary 
constitutive  assumption  in  damage  modeling  (Krajcinovic,  1996).  Further,  J-  is 
expressed  in  the  following  simple  form  (Lubarda  and  Krajcinovic,  1995) 

=  +  (64) 

where  Fi  and  F2  are  constants. 

Considering  the  case  of  isotropic  damage,  Pijki  can  be  taken  in  the  form  (49)  and 
(60)  should  reduce  to  the  simple  form 

iOab  =lF8ab  ■  (65) 


It  follows  that  Q  can  depend  only  on  i.e.  on  the  trace  of  fiab,  so  that 


dg 

dgiab 


Sab  ■ 


Consider  the  case  of  a  dilute  concentration  of  microcracks  where  the  interaction 
between  microcracks  as  well  as  the  effect  of  accumulated  microcracks  on  further 
microcrack  formation  can  be  neglected.  In  this  case,  the  rate  of  damage  evolution 
LJ  may  be  taken  to  be  independent  of  the  accumulated  damage  u>.  Dependence  of 
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u)  on  uj,  however,  is  implicit  in  Vt^c-  Thus,  to  eliminate  dependence  of  w  on  a;  the 
function  Q  is  taken  to  be  linear  in  ficc 

Q  =  (66) 

where  G  is  a  constant.  Expression  (66)  can  be  thought  of  as  the  first  term  in  a 
Taylor  series  expansion  of  Q  (the  constant  term  in  such  expansion  is  immaterial 
to  the  present  formulation).  Thus  addition  of  higher  order  terms  in  the  expansion 
introduces  dependence  of  a;  on  u;. 

Under  damage  loading  conditions  ci;  >  0  and  the  left  hand  side  of  the  dissipation 
inequality  (44)  becomes 

♦  2 

^  Ar7r  +  (1  io)  Cj  Cijkl  £ij  ^kl  ■) 

T 

which  is  always  positive  since  Cijki  is  positive  definite.  Thus  the  requirement  of 
positive  dissipation  is  identically  satisfied.  Substitution  of  (64)  and  (66)  into  (60) 
and  use  of  (65)  yield  damage  evolution  in  the  form 

lli  =  a  e'ij  e'i^  +  ^  Skk  £kk  ,  (67) 

where 

ct  =  2  jPj  G  and  =  2  F2  G 

are  free  parameters  that  need  to  be  determined  from  the  damage  evolution  pattern 
in  a  considered  problem.  The  first  term  on  the  right  hand  side  of  (67)  represents 
the  effect  of  the  deviatoric  part  of  the  behavior  on  damage  evolution,  wherea.s  the 
second  term  represents  the  effect  of  the  hydrostatic  part. 


6.  CONCLUDING  REMARKS 

In  this  article  a  thermodynamically  consistent  framework  was  proposed  to  model  the 
coupling  between  linear  viscoelastic  deformation  and  microcrack  damage.  The  effect 
of  damage  was  incorporated  into  the  constitutive  equations  in  a  form  consistent  with 
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the  well-known  effective  stress  and  effective  strain  concepts,  and  damage  evolution 
was  related  by  the  concept  of  damage  surfaces. 

Several  tensorial  ranks  for  the  damage  variable  can  be  employed  in  the  context  of 
the  present  formulation.  However,  for  simplicity,  in  this  article  damage  is  represented 
by  a  symmetric  second  rank  tensor.  This  representation  is  capable  of  simulating 
some  changes  in  material  symmetry  induced  by  microcrack  damage.  A  more  general 
representation  is  to  take  the  damage  variable  in  the  form  of  a  double  symmetric 
fourth  rank  tensor  tOabcd-  Such  a  representation  is  capable  of  simulating  general 
damage-induced  changes  in  material  symmetry  (Krajcinovic,  1996).  This,  however, 
complicates  construction  of  the  functional  forms  for  the  mapping  tensor  Pijki  and 
the  damage  potential  Q;  since  both  have  to  be  functions  of  double  symmetric  fourth 
rank  tensors  {i^abcd  in  the  case  of  Pijki  and  the  conjugate  thermodynamic  force  Ctabcd 
in  the  case  of  ^).  For  fourth  rank  tensors,  the  definite  forms  of  the  integrity  bases 
and  invariants  are  not  yet  well  established  (Zheng,  1994). 

An  important  remark  is  that  in  the  present  formulation  the  thermodynamic 
force  conjugate  to  damage  depends  on  the  viscoelastic  internal  state  variables  as 
can  be  seen  from  relations  (42)  and  (43).  This  differs  from  previous  formulations 
by  Schapery  (1981,  1996a, b),  where  the  thermodynamic  force  conjugate  to  damage 
is  taken  to  be  independent  of  the  viscoelastic  internal  state  variables  and  to  depend 
only  on  the  elastic  (instantaneous)  part  of  the  deformation.  The  argument  put  forth 
by  Schapery  is  that  based  on  a  viscoelastic  fracture  mechanics  analysis  (Schapery, 
1984)  it  was  found  that  the  driving  force  for  existing  cracks  is  independent  of  the 
viscoelastic  internal  state  variables.  However,  damage  evolution  occurs  not  only  by 
the  extension  of  existing  microcracks,  but  also  by  the  nucleation  of  new  microcracks. 
In  this  general  case,  it  is  expected  that  the  state  of  the  viscoelastic  deformation  in  the 
material  should  have  a  direct  effect  on  the  formation  of  new  microcracks.  Hence,  the 
thermodynamic  force  conjugate  to  damage  is  expected  to  depend  on  the  viscoelastic 
internal  state  variables.  This  is  also  in  agreement  with  the  internal  state  variable 
formulation  for  the  coupled  elastoplastic-damage  behavior  (e.g.  Ju,  1989;  Hansen 
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and  Schreyer,  1994). 
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